GOCE卫星重力梯度数据反演重力场的滤波器设计与比较分析

在地固系下，地球外部引力位的球谐级数展开式为：

式中，$ {G}_{\mathrm{G}\mathrm{M}} $为地心引力常数；$ R $为地球平均半径；$ n $、$ m $分别为展开的阶和次；$ {\stackrel{-}{C}}_{nm} $、$ {\stackrel{-}{S}}_{nm} $为规格化的引力位系数；$ r $、$ \theta $和$ \lambda $分别表示地心向径、余纬和经度；$ {\stackrel{-}{P}}_{nm}\mathrm{c}\mathrm{o}\mathrm{s}\theta $为规格化的缔合勒让德函数；$ N $为展开的最大阶数。

在局部指北坐标系（local north-oriented frame，LNOF）下，重力梯度分量计算公式如下[12]：

式中，$ x $、$ y $、$ z $为LNOF下的坐标；$ {V}_{xx} $表示引力位$ V $对$ x $的二阶偏导数，其他类似。

根据式（1）、式（2）可以建立重力梯度张量与引力位系数的线性关系，并采用空域最小二乘法[5-6]估计引力位系数参数。该方法可直接采用沿轨重力梯度数据求解，并不需要对观测数据进行归算与格网化处理，其解算过程较为严密。假设$ \boldsymbol{X}=\{\mathrm{\Delta }{\stackrel{-}{C}}_{nm};\mathrm{\Delta }{\stackrel{-}{S}}_{nm}\} $为未知的位系数参数改正向量，则LNOF下卫星重力梯度数据的误差方程为：

式中，$ {\boldsymbol{v}}_{\mathrm{L}\mathrm{N}\mathrm{O}\mathrm{F}} $为LNOF下的改正向量；$ {\boldsymbol{V}}_{\mathrm{L}\mathrm{N}\mathrm{O}\mathrm{F}}-{\boldsymbol{V}}_{\mathrm{L}\mathrm{N}\mathrm{O}\mathrm{F}}^{0} $为LNOF下的重力梯度观测值残差向量；$ \frac{\partial {\boldsymbol{V}}_{\mathrm{L}\mathrm{N}\mathrm{O}\mathrm{F}}}{\partial \boldsymbol{X}} $为设计矩阵。

由于GOCE重力梯度数据是在梯度仪坐标系（gradiometer reference frame，GRF）下观测得到的，并且$ {V}_{xy} $和$ {V}_{yz} $的精度比其他分量约低2个量级，若将重力梯度张量从GRF转换至LNOF，则这两个低精度的梯度分量将会导致其他分量的精度降低。为了避免坐标转换引起的精度损失，可将式（3）转换至GRF，从而得到：

式中，$ {\boldsymbol{v}}_{\mathrm{G}\mathrm{R}\mathrm{F}} $为GRF下的改正向量；$ {\boldsymbol{V}}_{\mathrm{G}\mathrm{R}\mathrm{F}}-{\boldsymbol{V}}_{\mathrm{G}\mathrm{R}\mathrm{F}}^{0} $为GRF下的重力梯度观测值残差向量；$ \frac{\partial {\boldsymbol{V}}_{\mathrm{G}\mathrm{R}\mathrm{F}}}{\partial \boldsymbol{X}} $为设计矩阵，计算如下：

式中，$ {\boldsymbol{R}}_{\mathrm{G}\mathrm{R}\mathrm{F}}^{\mathrm{L}\mathrm{N}\mathrm{O}\mathrm{F}}={\boldsymbol{R}}_{\mathrm{G}\mathrm{R}\mathrm{F}}^{\mathrm{I}\mathrm{R}\mathrm{F}}{\boldsymbol{R}}_{\mathrm{I}\mathrm{R}\mathrm{F}}^{\mathrm{E}\mathrm{R}\mathrm{F}}{\boldsymbol{R}}_{\mathrm{E}\mathrm{R}\mathrm{F}}^{\mathrm{L}\mathrm{N}\mathrm{O}\mathrm{F}} $，表示LNOF至GRF的转换矩阵，其中，$ {\boldsymbol{R}}_{\mathrm{G}\mathrm{R}\mathrm{F}}^{\mathrm{I}\mathrm{R}\mathrm{F}} $为惯性基准坐标系（inertial reference frame，IRF）至GRF的转换矩阵，可由EGG_NOM_2数据提供的卫星姿态四元素计算；$ {\boldsymbol{R}}_{\mathrm{I}\mathrm{R}\mathrm{F}}^{\mathrm{E}\mathrm{R}\mathrm{F}} $为地固坐标系（earth-fixed reference frame，ERF）至IRF的转换矩阵，可从SST_PSO_2I数据中获取；$ {\boldsymbol{R}}_{\mathrm{E}\mathrm{R}\mathrm{F}}^{\mathrm{L}\mathrm{N}\mathrm{O}\mathrm{F}} $为LNOF至ERF的转换矩阵，具体形式为[12]：

重力梯度观测数据包含有色噪声，因此需要设计白化滤波器对式（4）的观测值残差向量和设计矩阵的各列进行卷积，而滤波处理后的式（4）就可以采用最小二乘法求解位系数参数。由于GOCE数据的极空白问题及重力场信号向下延拓的不稳定性，其法方程的求解是一个严重的病态问题。采用Kaula正则化方法[6, 12]可以有效提高求解的稳定性，为确定最优的正则化参数，以及考虑到$ {V}_{xx} $、$ {V}_{yy} $、$ {V}_{zz} $和$ {V}_{xz} $分量滤波处理后给定的精度不一定准确，本文采用方差分量估计方法[13]进行迭代求解，其法方程为：

式中，$ i $表示不同的重力梯度分量；$ {\boldsymbol{A}}_{i} $、$ {\boldsymbol{P}}_{i} $和$ {\boldsymbol{l}}_{i} $分别为第$ i $个梯度分量对应的设计矩阵、观测值权阵和残差向量；$ \boldsymbol{K} $为Kaula正则化矩阵；$ \widehat{\boldsymbol{X}} $为位系数改正向量的估值；$ {\sigma }_{i}^{2} $为第$ i $个梯度分量的方差；$ {\sigma }_{K}^{2} $为正则化参数对应的方差。

文献[7]指出重力梯度数据误差具有有色噪声特性以及存在与轨道周期有关的系统误差，并且在测量带宽内重力梯度$ {V}_{xx} $、$ {V}_{yy} $分量的误差约为10 mE，$ {V}_{zz} $、$ {V}_{xz} $分量的误差约为20 mE。以2009-11-02的GOCE重力梯度观测数据为例，图 1给出了$ {V}_{zz} $分量的误差PSD。在测量带宽内（图 1两条黑色虚线之间）误差约为20 mE，并且在低频1 CPR（约1/5 400 s）及其倍数频率处存在明显的系统偏差。因此，利用GOCE重力梯度数据反演重力场除了要设计合适的滤波器对有色噪声进行处理外，还要考虑重力梯度数据低频系统误差的影响。

图  1  GOCE卫星重力梯度$ {V}_{zz} $分量的误差PSD

Figure 1.  Error PSD of Gravity Gradient Component $ {V}_{zz} $ of GOCE Satellite

通过MA方法可以获取时间序列中的低频部分。因此，将重力梯度误差时间序列减去其移动平均值就能削弱低频系统误差。基于该方法，可以设计MA滤波器，它是一种等效的高通滤波器，对应的线性差分方程如下：

式中，$ c\left(n\right) $表示滤波处理后输出的误差时间序列；$ e\left(n\right) $表示输入的重力梯度误差时间序列；$ d $表示滤波窗口，本文根据测量带宽下界选取的滤波窗口为60 s。与通常的高通滤波器相比，MA滤波器具有计算量低和便于编程实现的特点。另外，为避免滤波可能带来相位漂移的影响，该滤波器要同时处理式（4）的观测值残差向量和设计矩阵。

低频系统误差主要出现在1 CPR及其倍数频率处，因此可考虑获取重力梯度误差时间序列在这些频率的拟合值，该拟合值可看作低频系统误差，再将其从原始误差时间序列中减去，从而削弱低频系统误差的影响。基于该方法可设计CPR经验参数滤波器，其拟合低频系统误差的方程为[14]：

式中，$ r\left(n\right) $为拟合的低频系统误差时间序列；$ \omega \approx 2\mathrm{\pi }/5\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }400 $ s，为卫星轨道角速度；$ {x}_{0}, {x}_{1}, \cdots , {x}_{2i+1} $为经验参数，可通过最小二乘在每个轨道周期内估计一组，$ k $最大取27（即频率为27 CPR，对应测量带宽的下界）。需要指出的是，在实际滤波处理时要对$ k $的取值进行进一步分析，这是由于式（9）在吸收低频系统误差的同时也会削弱重力场信号。另外，设计矩阵中没有低频系统误差，所以只需要对观测值残差向量进行滤波处理。

以$ {V}_{zz} $分量为例，图 2给出了MA滤波器处理后的误差PSD。可以看出，滤波处理后测量带宽内及高频部分的误差保持不变，而低频部分的误差被有效削弱。同时，该滤波器只是将低频误差削弱至测量带宽的误差水平，并没有将其完全滤除，这样可以避免滤除低频重力场信号。

图  2  MA滤波器处理后的$ {V}_{zz} $分量误差PSD

Figure 2.  Error PSD of $ {V}_{zz} $ Processed by MA Filter

同样以$ {V}_{zz} $分量为例，图 3给出了CPR滤波器中参数$ k $分别取2、5、8（2 CPR、5 CPR、8 CPR）时，滤波处理后的误差PSD。可以看出，在对应CPR频率处，低频误差被有效削弱，但也削弱了CPR频率之间的误差，并且随着$ k $的增大，低频部分削弱更为明显。因此，为避免完全滤除低频部分，$ k $值不宜过大。此外，当$ k $取5和8时，高频部分保持不变，而$ k $取2时，高频部分出现明显的误差，并且当$ k $值过小时，低频系统误差的滤除效果不足。综合以上分析，后续$ k $值均取5。

图  3  CPR滤波器处理后的$ {V}_{zz} $误差PSD

Figure 3.  Error PSD of $ {V}_{zz} $ Processed by CPR Filter

针对重力梯度误差中的有色噪声，本文基于ARMA模型构建了ARMA去相关滤波器，其线性差分方程如下[15]：

式中，$ c\left(n\right) $表示输入的有色噪声序列；$ w\left(n\right) $表示ARMA滤波器处理后输出的白噪声序列；$ p $和$ q $为ARMA模型的阶数，可通过赤池信息准则（Akaike information criterion，AIC）来确定；$ {a}_{i} $和$ {b}_{j} $为滤波系数，可通过ARMA模型拟合$ c\left(n\right) $来确定。与MA滤波器类似，ARMA滤波器需要同时处理观测值残差向量和设计矩阵。

图 4为MA+ARMA级联滤波处理后的Vxx、$ {V}_{yy} $、$ {V}_{zz} $和$ {V}_{xz} $重力梯度分量误差PSD。滤波后4个分量的误差均表现出明显的白噪声特征，则此时式（7）中的权阵可近似为对角阵直接进行求解。

图  4  MA+ARMA级联滤波器处理后4个重力梯度各分量的误差PSD

Figure 4.  Error PSD of 4 Gravity Gradient Components Processed by MA+ARMA Cascaded Filter

图 5为CPR+ARMA级联滤波处理后的$ {V}_{xx} $、$ {V}_{yy} $、$ {V}_{zz} $、$ {V}_{xz} $重力梯度分量误差PSD，滤波后4个分量的误差同样表现出明显的白噪声特征，但其低频部分与图 4有差异，这是由于CPR经验参数滤波器滤掉了更多低频系统误差。

图  5  CPR+ARMA级联滤波器处理后4个重力梯度分量的误差PSD

Figure 5.  Error PSD of 4 Gravity Gradient Components Processed by CPR+ARMA Cascaded Filter

由于GOCE数据和GRACE数据在重力场反演中具有频谱互补性[12]，同时为了与已发布的联合模型进行比较，本文将联合GOCE数据和GRACE数据反演重力场模型。将这两类数据视为不同类型的观测值，可采用最小二乘联合平差方法估计引力位系数，其实质是将GOCE数据和GRACE数据形成的法方程进行叠加求解。式（7）为GOCE数据法方程，由于GOCO01S的研制采用了ITG-GRACE2010S[16]的法方程进行联合求解，为便于比较，本文也采用ITG-GRACE2010S模型的法方程进行联合求解。此外，在法方程叠加的过程中要合理确定GOCE数据和GRACE数据的权重，由于无法获取GRACE数据的观测方程（不能采用方差分量估计进行定权），本文采用经验定权的方法。